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. Abstract. The exploration of cold polar molecules in different geometries is a rapidly 

^ ' developing experimental and theoretical pursuit. Recently, the implementation of 

optical lattices has enabled confinement in stacks of planes, the number of which is 
also controllable. Here we consider the bound state structure of two polar molecules 
' confined in two adjacent planes as function of the polarization angle of the dipole 

moment of the molecules. We prove analytically and present numerical evidence for the 
existence of bound states for arbitrary dipole moments and polarization directions in 
' this two-dimensional geometry. The spatial structure of the bound states is dominated 

O . by two-dimensional s- and p-waves, where the latter exceeds 40 percent over a large 

range of polarization angles for intermediate or strong dipole strength. Finally, we 
consider the influence of the dimer bound states on the potential many-body ground- 
state of the system. 
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1. Introduction 

A strong experimental drive in the field of polar atoms and molecules has realized 
controllable samples in the rotational and vibrational ground-state that are close to 
quantum degeneracy [H El El Hj |5l E]. These heteronuclear systems have a number 
of very interesting properties due to the long-range and anisotropic dipole-dipole force 
which can give rise to highly non-trivial many-body states in both the weak- and strong- 
coupling regime [Tj [8] . The attractive head-to-tail configuration can, however, lead to 
strong chemical reactions [6] or many-body collapse of the system [9], and confinement 
in optical lattices has been suggested as a means of avoiding this problem ^0\. These 
confined one- or two-dimensional geometries have led to a number of predictions of 
novel few- and many-body states [T0l[IIl[T2l[T3l[Tl[T5l[T6l[I71[Ii[l9l|2^ 
very recently the first experimental implementation of a multilayered stack of pancakes 
containing fermionic polar molecules was reported [23] . 

Here we consider the case of two adjacent layers. However, even in this seemingly 
simple case there is a competition of intra- and interlayer interactions which can vary 
between repulsion and attraction as one changes the polarization angle of the dipole 
moments with respect to the layers. In the present paper we will be concerned with 
few-body states with one particle in each layer in order to describe the simplest complex 
in such a system in detail. The case of dipoles oriented perpendicular to the layers was 
considered from the few-body bound state and scattering point of view in previous works 
[211 EHl EHl ET] . At the so-called 'magic' angle where the intralayer repulsion vanishes in a 
one-dimensional trap the few-body bound state structure was also discussed [28 1 [29 | [30] . 

To our knowledge, the full two-body bound-state problem as a function of the 
polarization angle and the dipole moment has not been studied previously. This problem 
is highly non-trivial due to (i) the anisotropy and (ii) the vanishing integral over space of 
the potential for arbitrary polarization angle. The problem is a specification of a more 
general problem of two particles in two dimensions interacting via anisotropic potentials 
where the net volume is negative or zero. The more general problem has been addressed 
only in a recent letter [31], where the details of analytical and numerical methods are 
not elaborated. 

In this paper, we present in section 2 analytic results for the dipole-dipole potential. 
They are specialization of the derivation in Appendix A for an arbitrary potential. In 
section 3, we describe the novel numerical method based on stochastic variation along 
with the computed results for the dipole-dipole potential. We present energies, wave 
functions, and expectation values of relevant operators as functions of the potential 
strength and the polarization angle. One of our main results, analytical as well as 
numerical, is that the bilayer system has a bound state for any polarization angle and 
any value of the dipole moment. We also calculate a partial- wave decomposition that 
characterizes the geometric structure of the wave function which indicates the likely 
symmetries of the corresponding many-body problem. In section 4 we present a first 
application of our results in a many-body context. We consider the limit of strong 
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coupling where the system forms bound bosonic dimers that can potentially form a 
(quasi)-condensate. Finally in section 5 we brifly summarize and conclude. 

2. Analytic results 

The general setup we consider consists of two particles confined to two spatial dimensions 
with a pair potential V{r) depending on the relative coordinate f. In polar coordinates, 
r = {x,y) = {r cos if , r simp) , we have the Schrodinger equation: 



where is the wave function, the reduced mass, A is the dimensionless strength, 
AV^ = 2jj,d^V/lT? and = 2fj,d'^E/h'^, V is the potential and E the energy, d the unit of 
length, and s = r/d, is the reduced coordinate. 

The investigation of possible bound states was briefly sketched analytically in 
[3T] for an arbitrary potential in two spatial dimensions (2D). The result is that any 
cylindrical or non-cylindrical potential has a least one bound state provided the volume 
of the potential is negative or zero. We give more details of a similar derivation in 
Appendix A. The difficulties are centered around exceedingly weak potentials where 
no ordinary perturbation treatment is applicable because no unperturbed solution 
is available. The binding energy approaches zero as the potential vanishes and the 
continuum is approached. Thus the limiting energy is zero and the corresponding wave 
function is uniformly distributed over all space. 

We specialize to the dipole-dipole potential arising for the system of two polarized 
molecules of mass M confined to two parallel planes separated by a distance d as 
shown in figure [H The corresponding dipole-dipole potential, V, projected to this 
two-dimensional geometry is 



where D is the dipole moment and 6 denotes the polarization angle measured from 
the layer plane to the z-axis which intersects the two layers at right angles. 

The potential in equation ([2]) is found in the ideal limit of zero-width layers. 
The dipole polarization is measured such that for 9 = 7i/2, the dipoles are oriented 
perpendicular to the layers as in [211 [26l [27]. One can take corrections to the zero- 
width layer limit into account by integrating out a gaussian in the transverse direction. 
However, the corrections are second-order in the width, w, and we neglect them as we 
here are interested in the w ^ d limit. 

We have to solve the 2D Schrodinger equation in equation ([1]) with the potential 
in equation ([2|). The reduced mass is /i = M/2 and the dimensionless dipolar strength 
is given hy U = MD"^ /{fi^d), which is a measure of the ratio of potential to kinetic 

I In SI units we have = (P /Aireo when d is dipole moment of the molecules. However, in this paper 
we use d to denote the interlayer distance. 



1 d d 1 



(1) 



s ds ds dip 






Figure 1. Illustration of the setup consisting of two dipolar particles of mass M 
moving in parallel planes separated by a distance d. Their dipole moments, D, are 
assumed to be aligned by an external field at angle 6 with respect to the planes. 



energy. We will also consider the case where [/ < which is also physically realizable as 
explained below. In the notation above we have X = U and we will use these notations 
interchangeably in order to emphasize the generality of the analytic approach presented 
in Appendix A. 

This potential is invariant under reflection in the x-axis and has the peculiar 
property that / dxdyV{x,y) = for any 9. In particular, it does not fulfll the 
Landau criterion for bound states in two dimensions [32], which states that a bound 
state always exists for J dxdyV{x,y) < 0. An early existance proof was given in [33] 
using a method that is not well-suited for expansions in the strengh of the potential. A 
discussion of such an expansion appeared in [27] but only for case where 9 = 7i/2 and 
cylindrical symmetry holds. Here we are interested in the appearance and properties of 
bound states for arbitrary 9. A partial-wave decomposition of the potential in the basis 
{1, cos cos(2y9)} (which are the only non-zero terms) leads to 

XV{s,(p) = Vo{s) + Vi{s)cos(p + V2{s)cos{2(p) (3) 
N s^cos^^ 

which we will refer to as monopole, dipole, and quadrupole terms, respectively. The 
monopole potential Vq, has in itself zero net volume and it vanishes identically for 9 = 9^ 
where sin^^c = 1/3- The dipole term only vanishes for = and 7r/2, whereas the 
quadrupole term is flnite except at = 7r/2. Thus for 9 > 9^ and U > 0, the monopole 
term has an inner attractive pocket and a repulsive barrier outside s = \/2, and vice 
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versa for 6 < Oc- For [/ < the story is reversed. We expect the monopole term to be 
most important for the system properties, at least when it is non-vanishing away from 
6 = 6c- However, the monopole term is, except for the factor of (3 sin^ ^ — 1)/2, identical 
to the full potential at 7r/2, i.e. we know from previous work that it always supports 
bound states [261 EH [33] . We also know that the configuration with an attractive inner 
pocket and a repulsive outer barrier leads to considerably stronger binding than in the 
reversed case [26]. We will see this explicitly in the energies presented below. 

It is very important to notice that the angle 9c is different from the magic angle, 
6'*, where the potential of two dipoles moving in one dimension vanishes (determined 
by cos^ 9* = 1/3) [8j. This demonstrates an important difference between one- and two- 
dimensional dipolar systems. We will address this fact in more detail when we discuss 
many-body physics below. 

First we specialize the analytical results in Appendix A derived for general 
interactions to the dipolar potential. A partial wave expansion is employed in analogy 
to that of the decomposition in equation ([3]). For the dipole potential in equation ([2]) 
the resulting energy expression from equation (IA.25|) then becomes 
T, / 2{1 + UB, + U'B2) \ 

where 7 is Euler's constant and the coefficients, v4o,v4i and Bi,B2, are defined by 

Ao =^M'^+^sm'29 + ^cos'9, (8) 

Ai = + 0.0053 sin^ 2^ cos^e - 0.0033 sin 2^ cos^^ 

- 0.0019 cos^ 9 - Mc(0.0349 sm^{29) (9) 
+ 0.0054 cos^(e) + 0.0156McCos2 9 + 0.0343M^2) , 

5i = - 1.204M, - — cos^ 9 , (10) 
16 

B2 = 0.8382Mc{Mc + 0.0667 cos^ 9) (11) 

- 0.0037 sin^ (2^) + 0.0894 cos^^ , 

Mc = lsin\9)-^-. (12) 

The expression for A2 is much more elaborate consisting of more than a hundred terms 
each given as double, triple, and quadrupole integrals over well defined functions. We 
refrain from showing them all here. Note that since the spatial integral of the potential 
vanished, the term A2 has to be considered in an expansion to second order in U (see 
also 

The only non-zero matrix elements for the dipole potential in equation ([2]) are 
Vm,m±2- This implies that the m = 0, 1, and 2 partial waves are sufficient 
to get all contributions to the orders given on the right hand side of equation ([7]). Higher 
partial waves beyond m = 2 only contribute to the wave function for the dipole potential 
through orders of U that are higher than those in equation [71 

The scattering length, a, is usually considered to be a measure of the most crucial 
model independent property of any potential. In the present case it is a function of 
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strength and polarization angle, and to any order related to the energy as in [31], i.e. 
through the general equation 

^ = -2— exp(-27) . (13) 

Any accuracy of E is then directly transferable to a through equation ( IT3l) . The energy 
can be calculated to any order in powers of U, as second order in equation ([7]). Then 
the scattering length becomes 

where the dependence on strength and polarization angle now is explicit. 

The energy, and scattering length, very close to threshold is exponential in 
as seen in eqs.© and ([HD, and determined by the polarization angle through Aq. The 
first order terms, {Ai, Bi), in U exhibits the difference in approach to threshold for the 
different signs of the strength, U. The second order terms, {A2,B2), are necessary to 
get the correct ^/-independent pre-exponential factor in the energy. Here B2 is both 
much simpler and more significant than A2 which consists of sums of a large number of 
contributions expressed as definite integrals. 

The expressions simplify substantially for 6* = 7r/2. In ^35i| the energy is calculated 
for 6* = 7r/2 to the order given in equation ([7]), including the A2 term, in agreement with 
our result. The computation of the energy in |27j for 9 = tt/2 deviates from our results 
in equation ([7]) in the first order correction. 



3. Numerical procedure 

The potential is in general anisotropic and the wave equation is not easy to solve by 
discretization or integration. We therefore turn to the stochastic variational approach 
using gaussian wave functions which has been successfully applied to other interactions 
[36] . However, in the limit of weak binding the wave functions become very small and 
spatially extended without structure at large distances. The special method to achieve 
convergence with a fair amount of gaussians is described in this section, and the results 
for energies and wave functions presented in the next section. 



3.1. Method for weakly hound systems 

For numerical calculations we employ the correlated Gaussian method which has been 
sucessfully used in a range of few-body problems in atomic physics [371 EHl EH] . The 
wave function \E'(x, y) is found through the variational principle as a linear combination 
of basis-functions Gj(x, y), 

-^basis 

^{x,y) = CiGi{x,y), (15) 
1=1 
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where A'basis is the size of the basis, q are the hnear variational parameters, and the 
basis-functions are chosen in the form of shifted correlated Gaussians, 

Gi{x,y) = e-^''-'^^^^^^''-'^\ (16) 

where the superscript T denotes transposition, q = (x, y)'^ is the column-vector of the 
coordinates, and where the elements of the symmetric positive correlation matrices Ai 
and the shift vectors Sj are the non-linear variational parameters. The explicit shifts 
employed here enhance greatly the flexibility of the correlated Gaussians specifically for 
the non-rotationally symmetric system at hand. 

According to the variational principle the wave function is found by minimizing the 
expectation value of the Hamiltonian, 

= \' (17) 

The linear parameters q are determined by solving the generalized eigenvalue problem 
He = E[^]^c , (18) 

where Ti and M are A^'basis x Abasis matrices, 

-Hij = {G^\H\G,) and = , (19) 

where c is the column- vector of the linear parameters c,. The overlaps and the matrix 
elements of the kinetic energy operator in equation (JTOll are calculated using the 
analytical expressions in Appendix B. The matrix elements of the interaction potential 
are calculated numerically using adaptive Gauss-Kronrod quadratures [ID]. 

The generalized eigenvalue problem equation (ITSl) is solved by first performing 
Cholesky decomposition of the (symmetric and positive definite) matrix JV = LL^ , 
and then solving the ordinary symmetric eigenvalue problem. 

He = E[^]d , (20) 

where H = L^^HL^^^ and c = L^c, using standard linear algebra methods |40j . 
The non-linear parameters - the elements of the correlation matrices Ai and the shift 
vectors si - are determined through stochastic sampling. Basically, the minimum of 
the energy is found by sampling a large number of random sets of Abasis Gaussians 
with randomly generated correlation matrices and shift vectors. The elements of the 
correlation matrices are generated as ±1/6^ and the elements of the shift vectors as ±6, 
where 6 is a stochastic variable with the dimension of length sampled from an exponential 
distribution with certain length parameter /. To ensure the positive definiteness of the 
correlation matrix the following procedure is used: First a diagonal matrix with positive 
elements 1/6^ is generated, and then a random orthogonal transformation is performed. 
For ordinary quantum systems the length parameter, /, of the exponential distribution 
should be chosen close to the typical range, d, of the interaction potential between 
particles. For the dipolar potential considered in this paper, it is in fact the interlayer 
distance d which explains the abuse of notation. If the binding energy B of the system 



is of natural size, that is ~ such sampling proves very effective. However, for 
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weakly bound systems, where B <^ ^^pi two different length scales are of importance: 
One is the interaction range d and the other is the typical size |«|~^ of the tail of the 

wave function (where a = isj^^)- We therefore subdivide the basis into Ai'short, short- 
range, and Mong, long-range, Gaussians, iVbasis = ^short + A^iong, where the short-range 
Gaussians are sampled from the exponential distribution with range d and the long- 
range Gaussians are sampled from the exponential distribution with range The 
latter is calculated from the current estimate of the binding energy. 




-0.02 

5 10 15 20 25 30 35 40 45 

Figure 2. Illustration of the convergence of the energy with respect to the basis size: 
the energy E of the ground state as function of the size A'^basis of the basis. 



Since the long-range Gaussians are introduced specifically to better describe the 
asymptotics of the wave function, they can be chosen in a much simpler form, 

G = e"^'/''' , r = ^x^ + 1/2 . (21) 

The typical convergence plot of the binding energy as function of the basis size is 
presented in figure O It shows that the energy is converged to within four significant 
digits with the basis size of about 45 Gaussians. The convergence of the tail of the 
wave function is illustrated in figure |3] where it is compared to the analytic asymptotic 
form. Clearly, addition of the long-range Gaussians significantly improves the quality 
of the wave function at very large distances. In this example about 46 Gaussians, 30 
short-ranged plus 16 long-ranged, can accurately describe the the wave function up to 
the radius of about SOOrf. 



4. Results 



The general derivation in Appendix A demonstrates that there always is at least one 
bound state in 2D for any anisotropic potential with zero net volume as obtained in 
a different manner in [33]. When U becomes small we expect universal behavior of 



Bound dimers in bilayers of cold polar molecules 

1 



0.1 



0.01 



\const X -^"0(10; 

30 + 
30 + 12 
30 + 8 



U = 0.9; e = 9* 




40 600 800 

(v/^^+F)/t^ 



1000 



Figure 3. Illustration of the convergence of the radial wave function with respect 
to the basis size: the radial function ^o{r) from equation (jA.17|) as function of s = 
\/ + y^/d calculated with bases of different sizes: iVshort = 30 and iViong = 8, 12, 16. 
For comparison the analytic form of the tail, -Kodalr), is also shown. 



energies and radii [3T| IMf HI] . Using the stochastic variational approach in the small 
U limit, our results for = 7r/2 approach the universal behavior of the energy which to 
leading order scales like \n{—E) oc —l/U"^ as discussed previously [261 123 1211 [33]. For 
other values of 6 we expect the same scaling for very small U, however, the range of U 
around zero where this applies has a strong dependence on 9 as seen by the energies 
presented below. We also expect differences for general 9 between positive and negative 
U in the limit ?7 — )■ as for ^ = 7r/2 |26]. The question of how the binding energy 
approaches universality is investigated in more details in [31] . 

The energies have been calculated using the correlated gaussian approach. In figure 
Hjwe exhibit the results as a function of f/ > for a selection of polarization angles. 
At small U the energy decreases very fast with decreasing U as noted already in [33] , 
whereas at larger U we find a linear dependence on U as argued in [26] for 9 = 7c/2. 
The binding energies decrease dramatically as 9 approaches zero. We stress that we 
numerically find a bound state for any value of U also in the particularly unfavorable case 
of ^ = as well as for 9c (sin^ 6'^ = 1/3) where the decisive monopole potential vanishes 
identically. Notice that there are more bound states for larger U but we restricted our 
discussion to the single bound state regime. 

The U < case is also of great interest as that potential can be generated by using 
microwave-dressed molecules. In [151 [El [22| an AC light field directed perpendicular 
to the layers was used to create the 9 = 11/2 potential with U < 0. A straightforward 
calculation shows that if the laser hits the layers at an angle, 9, the potential is the same 
as for a homogeneous electric field at angle 9 but with negative U. For U < we again 
find numerically that for all values of the strength the two-body system has bound states. 
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Figure 4. Bound state energies as a function of U for different angles. Insets show 
contour plots of the potentials with valleys in bright (blue) and hills in dark (red) 
colors. 

The results for the binding energy at different angles are shown in figure [5] as function of 
\U\. The first thing one notices is that the overall magnitude of the bound state energy 
is smaller than that for U > 0. At 6 = ir /2 this can be understood as the potential 
has a repulsive core at = 0, forcing the state to reside in the shallow attractive 
pocket at intermediate distance. In turn, this gives a much smaller binding energy. 
This qualitative behavior of the potential persists until 9 decreases below 9c where the 
monopole changes sign. Then the potential changes overall character to become more 
attractive with inner attractive pocket and outer repulsive tail. The U < results thus 
show maximum binding at ^ = which is, however, still about a factor of three smaller 
than the U > case at its most favorable angle of 9 = tt/2. 

Numerical and analytical results are compared in figure [61 To get the most accurate 
comparison we show the exponent in equation multiplied by the square of the 
strength, U'^. The approach in the figure of analytical and numerical results is reassuring 
in the limit of relatively weak potentials. Since we neglected the y42-contribution this 
approach indicates its minor significance. 
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Figure 5. Same as figure H] but for U < 0. Note tlrat tlie vertical scale is different 
from that of figure [H 

4.I. Wave functions 

The structure of the bound state wave functions can be seen from the partial wave 
decomposition. The results are shown in figure [7] for a strong coupling of f/ = 10 and 
a weaker one of t/ = 4. The probabilities are normalized so that they sum to one. We 
note that the contribution of m > 2 is only a few percent with a maximum at 6* = of 
5% in m > 2 terms. As expected we find that m = becomes dominant for 6 — ^ 7r/2 
as we approach cylindrical symmetry. Interestingly, close to = we also find a very 
large m = component, no m = 1 content, and a significant m = 2 contribution. The 
remaining content of the wave function is found in the higher m contributions. The 
fact that m = 1 has no weight for ^ = can be understood from the symmetry of the 
potential. For x — —x the m = 1 term changes sign, whereas the potential is invariant. 
Interestingly, as we go away from 6 = 0, the m = 1 component raises rapidly and stays 
on the order of 40% until we reach 6 = O^. a.t which it starts to decline as for m = 2, in 
line with the restoration of cylindrical symmetry at 6 = tt/2. For U > 10 the m = 
component can be even more suppressed in comparison to m > for intermediate 6, 
whereas for positive U < A the m = component will eventually dominate as one 
approaches the universal limit discussed above. 
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Figure 6. The function F{U) = -U^ln{\E\Md'^ /{2h'^)) for different polarization 
angles 6. The dashed red curves are calculated numerically and the solid blue curves 
are from equation ([T]). The relatively small contribution from the yl2-term is neglected 
in this comparison. 



We have found similar results for the U < when taking into account that the 
angle 9 for U > correspond to angle — 9 ioi U <0 and vice versa. This is in fact 
an exact symmetry of the dipole part of the potential and an approximate one for the 
monopole term since 9c is close to 7r/4. For U = — 10 we find that there is a window 
6'c < < 1.1 in which the m = 1 term is around 40%. Interestingly, we find that the 
partial- wave content for f/ < is almost exclusively m = and m = 1. This is perhaps 
surprising as the potential in the m = 2 channel is non-vanishing except at ^ = 7r/2. 

The potential has completely different forms for different polarization angles as 
illustrated on the contour plots in the inset of figure HI For 9 = n/2, the potential 
is cylindrical while asymmetry appears for decreasing 6'-values, and eventually two 
minima emerge when 9 approaches zero. For small 9, the potential looks like a harmonic 



Bound dimers in bilayers of cold polar molecules 



13 



1 

0.9 
0.8 
0.7 
0.6 
0.5 

0.4 

0.5' 

0.4 



^1 1 1 


1 1 1 ■ ' 




^^^^ 

_ - •- 'J- 








^ - ■- — 




'''''^^^ TTi = 




1 1 1 1 



0.2 



0.4 



0.6 



0.8 



1.2 



1.4 



g 0.3 
CIh 0., 



0.2 0.4 0.6 0, 




1.4 



m = 2 



1.4 



1.6 



1 1 


1 1 1 1 


1 


" " " ~ - .^^^^ 171=1 


-1 X 




/ ✓ 
/ x 


^ ^^^^ 

•lb ^^^^ 


/ ✓ 

f ^ 


^ ^^^^ 


L «• 1 1 1 


1 1 1 



1.6 



1.6 



Figure 7. Partial waves probabilities, P,n — 7r(l + (5,„o) /q°° rfpp|^'m(p)P with 
^'(p, 0) = X^m *™(/°) '^*^^(''^'?^)' f'^'" t'^^ bound state wave function at J7 = 10 (solid 
black) and U = A (dashed red) as function of polarization angle 6* for to = 0, 1, and 2. 
The vertical line indicates 9 = On- Note the different vertical scales. 



oscillator along the y-axis for small x and y. The depth of this harmonic well around 
zero is about twice as large for ^ = and f/ < in comparison to the depth of the two 
wells in the x-direction for ^ = and f/ > that is shown in the inset of figure H] and 
figure [HI A simple gaussian wave function should therefore be a fair approximation to 
the full problem and in turn the lowest partial-waves should dominate. 

The wave functions for strongly bound states mimic the contours of the potential. 
As the strength, as well as binding energy, decreases the wave function is spreading 
out to larger distances and approaching cylindrical symmetry, as illustrated in figure [HI 
The probability decreases in all points of space and approaches zero uniformly outside 
the potential. However, inside the potential the shape of the potential is maintained 
even for vanishingly small strengths where the probability also approaches zero. This 
behavior is necessary to provide binding which arises from the attractive potential at 
small distances. In turn, the modified Bessel function, i^odals), is approached for 
vanishing strength which corresponds to a wavefunction that is roughly constant in 
space until the distance, s, is comparable to l/|a|. 
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e = Ti/A 





Figure 8. Contour diagram for different polarization angles of the dipole-dipole 
potential divided by U (upper part) followed by the probability distributions below. 
All lengths are in units of d. 
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5. The Many-Body Bilayer System 

The bilayer system has an interesting many-body structure with combination of 
attractive interactions that can induce pairing and repulsive interaction that tend to 
suppress such effects. This was discussed recently for the 6 = -n (2 case in [191 EO]- Here 
we consider the strongly-coupled limit (large U) where the bound two-body dimers are 
expected to be the relevant degrees of freedom. As the dimers are effectively bosons, they 
are capable of forming a (quasi)-condensate under the right conditions [SO] . However, as 
is well-known from BCS-BEC crossover studies [12], this is only expected to occur when 
the density is low. Unfortunately, the Berezinskii-Kosterlitz-Thouless (BKT) transition 
[ISj m] that governs this two-dimensional system has a critical temperature that is 
proportional to density [HI [20]. Therefore, a compromise where the dimer condensate 
occurs at not too low densities would be optimal to allow experimental access to this 
unusual many-body state. 

A further interesting complication is the fact that each of the layers can accomodate 
various coherent many-body states when we consider them independently. Proposals 
for the ground state include p-wave superfluids [HI [16] and density waves [17], [18]. 
Of course, when more than one layer is present the long-range inter- and intralayer 
interactions are competing, and at this point it is not clear what state is favoured for 
arbitrary directions of the polarization. It is known that the density-wave instability 
will be enhanced and occurs at smaller coupling strength for bi- and multilayer systems 
[2T] . We note that in the f/ < bilayer particle-hole coherence reminiscent of 

ferromagnetism has even been suggested [ilS]. Here we will be concerned mainly with 
the BCS-BEC crossover scenario in the strongly-coupled limit, but we will also estimate 
the appearance of a density-wave state and comment on possible superfluids. 

As the criterion for the onset of condensation of dimers, we consider the point at 
which the chemical potential becomes negative [20], i.e. 

^(t/, 6) = ^nV^siU, e) + Ep- ]^Eb{U, 9), (22) 

where Eb is the dimer binding energy and V^g is the long-wavelength (zero momentum) 
effective momentum-space interaction between two dimers. Here we include both the 
binding energy and the dimer-dimer interaction, and we also include a term for the Fermi 
energy, Ep, that the constituents of the dimer inherit from their layer. The density 
of dimers (equal to the single-layer density when the layers have an equal number of 
molecules) is denoted by n. To obtain the effective interaction, one must in principle 
integrate out the wave function of the dimer and include all inter- and intralayer two- 
body terms [20j. However, here we are only interested in the long- wavelength limit 
(momentum zero) in which the interlayer term vanishes [13]. This gives 

where ^2(2;) = (3x^ — l)/2. For the layer width, we take w/d = 0.2 in the following. 
Notice that V^g is attractive for 9 < 9c, vanishes at 9c, and repulsive for 9 > 9c- The 
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attraction for 6 < 6c results in a negative compressibility in a single layer P^. We stress 
again that 6c is much smaller than the angle at which the intralayer repulsion vanishes in 
a one-dimensional system. In this sense 6c is a special angle for the intralayer repulsion, 
whereas it has no dramatic effect on the binding energies which vary smoothly around 
6 = 6c- Combining the formula above, the final expression for /i becomes 

-/^ = in^ - P2{sin6) + 1 - —^Eb, (24) 



where we use the Fermi momentum kp = Ann for fermions in a single layer in place of 
n. 

The lines of /i = for selected angles are shown in figure |9] in the {U,kpd) plane 
for 1.5 < U < 5. For U > 2 the dimers have significant binding energy and can we 
treat them as localized bosonic objects. For 6 = tc /2 we present results both with and 
without the intralayer term which is clearly seen to shrink the region of potential dimer 
condensation. For 6 = 6c, the intralayer term vanishes and we find a larger region of 
/i < 0. For 6 < 6c the region would in principle become even larger, however, the 
intralayer term is then attractive and can lead to instability and collapse [I^. We 
therefore expect the line for 6' = 6'c to provide a boundary for how large the BEC region 
can become when tuning the angle within our approximations. In figure |9] we also show 
the line above which the density-wave instability appears at angle 6 = 6ciTL the random- 
phase approximation as discussed in Ref. [21] (above full blue line, denote by DW). To 
study the crossover outside the density-wave region, we see that low densities are indeed 
needed. 

We now consider the important question of the finite temperature behavior of 
the system. In the large U limit, the BKT transition temperature is maximal at 
kBTBKT = Ef/8 [151 [m EO], or 

Tb^t = 765 —^^nK. (25) 
10^ cm"^ M 

If we consider LiCs molecules (which can have a dipole moment of up to 5.5 Debye) 
with d = 0.5 /im, then the fi = phase for 6 = 6c at U = 10 is at = 8.1 ■ 10^ cm^, and 
thus Tbkt ~ 4.9 nK. While still very low, this is a significant increase over the sub-nano 
Kelvin temperatures for ^ = 71/2. 

The p-wave superfiuid for single layers is estimated to appear for angles sin 6 < 2/3 
[13]. However, density- waves occur for sin^ ^ > 1/3. This implies that there is a rather 
large regime (2/3 < sin 6* < 1) where the single-layer p-wave superfiuid should be absent 
but where BCS-BEC crossover induced by the interlayer interaction is possible. We also 
note that there is an intermediate region where both p-wave superfiuid and density- 
wave have been proposed as the ground state and one could imagine something like a 
supersolid phase as recently discussed in a two-dimensional optical lattice with fermionic 
polar molecules and perpendicular polarization [15|. In the strong-coupling limit we 
expect that the allowed few-body states will play an important role in determining the 
system properties, and calculation of the binding energies of states with three or more 
particles for arbitrary directions of the polarization are being pursued. 
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Figure 9. Lines of vanisliing cliemical potential for sin^ 6'c — 1/3 (solid black) and 
9 = tt/2 with (solid red) and without (dashed red) intralayer repulsion. Above the 
solid lines we expect condensation of dimers to occur (BEC), whereas in the lower 
right part a many-body paired BCS-like state (BCS) should be the ground-state of the 
system. Also shown is the line above which a density wave (DW) for 6* = 6'c is expected 
in the bilayer system (solid blue). 

We expect that the partial-wave analysis presented in the previous section can 
help indicate what symmetries are possible and relevant for the corresponding many- 
body problem in the large U limit. The problem is of course still that the intralayer 
term is attractive in the long-wavelength limit for 9 < 6c, and we thus expect that 
the most stable system require 9 > 9c where the decomposition of the wave function is 
entirely dominated by the m = term. However, close to 9c we still have a substantial 
171 = 1 contribution (immediately to the right of the vertical line in figure [7]). We 
therefore expect a region of interest in which an exotic many-body state with non- 
trivial symmetry like a p-wave dominated or mixed symmetry superfluid would emerge 
in the bilayer. These indications are consistent with a partial-wave decomposition of the 
potential itself [16]. Combining this information strongly suggests that there is a very 
interesting crossover from weak- to strong-coupling in the corresponding many-body 
system as recently discussed for the 9 = tt /2 case fi9\ [20] . Similar considerations hold 
for the U < case. 
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6. Summary and Outlook 

We have studied a bilayer system of dipolar molecules for arbitrary orientation of the 
dipoles with respect to the planes. The two-body bound state structure was calculated, 
including energies and partial-wave decomposition of the wave function as function of 
dipolar strength and polarization angle. We proved that there is always a bound two- 
body state in the system, irrespective of strength and polarization angle of the molecules, 
and also verified this fact numerically. We argued that this follows from the fact that for 
small strength, the wave function must reside outside the region where the potential is 
non-zero. The results apply irrespective of the sign of the interaction strength. Negative 
strengths invert the dependence of energy on the dipole angle such that perpendicular 
polarization angle has the smallest binding energy. The structure of the wave function 
is dominated by the monopole component which decreases with the strength of the 
interaction. Up to moderate strengths, the monopole component is always larger than 
50 percent while the dipole component accounts for most of the remaining probability. 

The conclusion that zero net volume potentials always have at least one bound 
state could perhaps be reached in other ways. First, approaching this limit from small 
negative net volume with one bound state strongly indicates that the bound state 
remains. Second, a perturbation argument is tempting, e.g. let us assume that the 
cylindrical monopole potential always has a bound state (as shown in [33]), and then 
treat dipole and quadrupole terms as perturbations. We thus extend from a Hilbert 
space entirely of s-waves to include also p and ci-waves. To second order in perturbation 
theory this always gives a negative contribution and hence more binding than for the 
monopole potential alone. However, closer scrutiny of such argements and their practical 
implementation reveal that in the limit of weak potentials the perturbations are always 
of the same order as the initial potential. This type of perturbative argumentation 
always fails. The deeper lying reason is that the energy is a non-analytical function of 
the strength for zero net volume potentials. This is seen in the expression for the energy 
where zero'th and first order in the strength contribute to the non-analytical structure 
at the continuum threshold. Only terms higher than second order in the strength are 
perturbations around this strong singularity which is a characteristic feature of zero 
volume potentials in two spatial dimensions. 

Implications for the many-body physics of a bilayer were discussed in the limit 
of strong-coupling where the two-body bound states are expected to be the important 
degrees of freedom. We conclude that the region where (quasi)-condensation of two- 
molecule dimers is likely to occur can be enhanced by tuning the angle of the dipoles. 
In particular, at the critical angle, sin^ 9c = 1/3, where the long-wavelength part of the 
intralayer interaction vanishes, we expect the conditions are most favorable for accessing 
this phase where dimers condense. This critical angle is different from the 'magic' angle 
at which dipoles moving on a line become non-interacting which has been discussed in a 
number of previous works [8], EH [29] . We also estimated the potential for a density wave 
instability in the bilayer [21], and demonstrated that BCS-BEC crossover and (quasi)- 
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condensation of dimers discussed here can occur also in an intermediate coupling and 
low density part of the phase diagram which is below the density wave regime. The 
possible roton instability in the bilayer system, as discussed for the perpendicular case 
in [20], is currently under active investigation. 

The results presented in this paper indicate that bound complexes of more than 
two particles must exist in the bilayer. In the case of one-dimensional tubes this was 
studied in some special cases in |28] and [29], while a more complete investigation of 
one-dimensional complexes as function of angles and dipole moment can be found in |30] . 
For the two-dimensional case, the method employed here can be extended to complexes 
of more particles and we plan such investigations in the near future. We also note that 
the external trapping potential that is present in each layer in experiments [23] can be 
easily accomodated in the current approach by introducing one-body harmonic oscillator 
terms. 

In conclusion, we find that bound states of dimers in a bilayer consisting of one 
particle in each layer are generic for particles interacting through the dipole-dipole force, 
irrespective of the dipole strength or polarization angle of the dipoles with respect to 
the layers. In general, the wave function contains several partial-wave components and 
therefore has interesting spatial structure. This suggests that few-body states with 
more than two particles will also have rich structure and it also implies that the many- 
body physics of the system is highly non-trivial. We sketched a phase diagram for the 
appearance of a dimer condensate as a function of polarization angle and showed that it 
is enhanced around the so-called magic angle. At this point the dimer contains a large 
admixture of higher partial waves and we expect the collective behavior of the system to 
reflect this fact. The many-body problem of a bilayer with polar molecules of arbitrary 
polarization angle therefore deserve further investigation. 
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Appendix 

Appendix A. General derivation 

We assume reflection symmetry in the for the potential in equation ([2]) but we 

shall otherwise proceed with the general derivation for any potential with this symmetry. 
The slightly more general formulation without any symmetry can be found in [31]. We 
decompose the wave function into partial waves, cos{mip), i.e. 




m=0 
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icos(m0)^(s,0), (A. 2) 



am.^m.(s) 1 



^/s (1 + 6mo)TT Jq 

where the corresponding contribution of the sin(m0) terms is zero due to the assumed 
symmetry and the coefficients are real. The normahzation of the functions $rn(s) will 
be chosen later. The rather artificial separation into coefficients and functions allows 
a simple procedure to compute higher orders in the strength A. The 2D radial wave 
functions, ^m{s), satisfy the system of coupled radial equations 

K + ^—J^^m + «'<J>^ = A V ^<!>iVUs) , (A.3) 
4s^ ^ am 

where the matrix elements, Vmi, carrying all information about the potential, are 

1 Z"^'" 

^rni = TT—, — ^ / cos{myD)cos{lip)V{s,ip)dip . (A.4) 

(1 + dmOjTT Jq 

For cylindrical potentials Vmi oc 6mi and the different m-values decouple. The regular 
radial solution to equation flA.31) at the origin provides the usual boundary condition for 
a centrifugal barrier potential, i.e. we choose the normalization of ^m{(^,s) such that 
lims_s.o s^^/^^'"<l'm(«, s) = 1. We inserted here exphcitly the dependence on the energy 
parameter, a, in We assume here the potential, and consequently also the right 
hand side of equation flA.3l) . diverges slower than when s — > 0. 

We use the Green's function formalism which in 2D is described for a cylindrical 
potential in ^7\. For anisotropic potentials we have more generally that the m- 
components, Om^m, are given by 

am^rn{a, s) = a^<l>„o(a, s) (A. 5) 

-AVaj / gm{a,s,s')Vmiis')^i{a,s')ds' . (A.6) 

The last terms vanish for s = and the boundary condition at s = is obeyed through 
the ffist term: 

= ^/^Jm{as){2/a^m\ , (A.7) 

where Jm is the Bessel function. The Green's function in the last terms of equation 
(1A.6I1 is expressed as 

gm{oi, s, s') = (A. 8) 

in terms of Hankel functions, Hm^ |1S]. At large distance equation (IA.6p . with the 
Green's function from equation (lA.Sp . may be rewritten as 

1 /2\™ 

am'l>m{(y,s) = -y/s i-j ml (A.9) 



X 



ml) 0,1 
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where the matrix elements, Cmi, axe: 

Cml = (A. 10) 

J 

Equation flA.Qp contains incoming and outgoing waves and is apphcable for scattering 
problems. For bound states, a has to be imaginary corresponding to negative energy, 
and Hm{(ys) in equation ( ]A.9I) diverges unless the coefficient vanishes, i.e. 



oo 



^Cmiai = -am- (A. 11) 

This is an eigenvalue equation, which only has non-trivial solutions for discrete values of 
a, and hence for the binding energy. However, the full radial wave functions, enter in 
the definitions of the matrix elements and the equations must be solved selfconsistently. 
We look for solutions in the limit of very weak strength, |A| <C 1, which implies that |a| 
(i.e. the energy) must also be very small. The m = component will then dominate in 
the solution to equation f lA.lip . because the centrifugal barrier suppresses higher partial 
waves. This suppression becomes more pronounced with decreasing binding energy. 
We now expand both coefficients, a^, and functions, $m, in powers of the strength, 

i.e. 

= Aa« + A2a(^) + ..., (A.12) 
$^ = $(^) + A$« + A2$(^) + ..., (A.13) 

where we leave the coefficient ao (without expansion) for normalization of the total wave 
function. In total we then find 



Oq 2m 

(2) 1 

ao 2m 



s'~"'Vrr,0is)ds , (A.14) 



POO 

/ s-^-"'V^ois)^'i\o,s)ds (A.15) 

— s'^-^V^,{s)^f\0, s)ds , (A.16) 
2m ^ ao Jo 



$i')(0, s) = -v^ /' s'Voois') \n{s'/s)ds' , (A.17) 

Jo ymo[s'){s'Y ™ds' 



^K.o(s')*i'Ho,.')-E4Ty^™(^')'^5°^(0'^'))d^'' 



Qim j^o 



Bound dimers in bilayers of cold polar molecules 22 

^l^\0,s) = - rds'goiO,s,s') 
Jo 

(1) 

X fVoois')^S\o,s') + y2^Vo.{s')<!>f\o,s')) , (A.20) 



^o(0,s,s') = v^ln^ , (A.21) 



gm(0, s, s') = — Vss'— — ^ . (A.22) 

The equations IA.17|A.19t and IA.20I provide the expansions for the coefficients and the 
wave functions in eqs. (IA.12P and equation (1A.13P when s ^ l/|ct|- The behaviour of 
the wave function at infinity is now given by the non-diverging piece in equation (IA.6p . 
i.e. 

hm s) ~ v^/j£)(as) , (A.23) 

or in the particular case of weak binding we have 

hm $^(a, s) ~ y/^H^\as)6^o . (A.24) 

s— >oo 

This behaviour is a consequence of the attractive and repusive centrifugal barriers for 

m = and m > 0, respectively. 

The energy can now be found to any order in powers of A. To second order the 

results can be written as 

2^2 

E= -_exp(-27) (A.25) 



X exp 



2(1 + ASi + A^Ss) 



-XI + A2(Ao + AAi + AM2) 
where the leading order constant is given by 



POO POO pZlT 

1=1 sVoQds= I I sV{s,^)dsd^. (A.26) 
Jo Jo Jo 

The potential dependent constants, {Ai, Bi), become increasingly more complicated with 
powers of A. When the net volume of the potential, XI < 0, is negative the weak binding 
limit for an arbitrary potential is given by the Landau expression [32] , which also turns 
out to be valid for anisotropic potentials with the appropriate definition of the volume. 
In the case of XI < 0, it is not necessary to retain A2 in (1A.25P to second order in A, 
whereas for / = it must be retained when calculating the corrections to the leading 
term. 

When / = corresponding to zero net volume of the potential the leading order 
term is given by Aq, i.e. 

POO 

Ao = - V^Voois)4'\o, s)ds + (A.27) 
^0 

poo 1— m roo 

T ^Vmois)ds v^Vo^(s')$L°Ho,s')ds', 

LnJo "'0 



m=/=0 
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where only first order from equation (1A.14P lias to be used in the derivation. 

This derivation is completely general for two-dimensional, anisotropic and reflection 
symmetric interactions. The symmetry requirement is only a minor simplification, 
and omitted in the derivation in [31]. The overall results are that there is always a 
bound state for very weak potentials with negative or zero net volume, and the weak 
binding threshold behavior of the energy is given by equation flA.25P to second order 



in the potential strength. The leading order term for zero net volume is given by 
equation ( 1A.27I) where only the first term contribute for cylindrical potentials since 
then Vom oc 5om- For non-cylindrical potentials even the leading order expression in 
equation ( 1A.27I) is rather complicated. 

Higher orders than those related to I and Aq are found for small energy by an 
iteration procedure through eqs.( 1A.6l) and ( lA.lOl) . The radial solutions are computed 



from equation ( 1A.6P which in turn are used to determine the C-matrix and the energy. 
This procedure can be repeated to give higher order corrections of both energy and 
wave function. Much care is necessary to include consistently all terms up to a given 
order because the resulting expressions contain many terms. The simplest is Bi which 
is found to be 



oo 



Bi= - Voo{s) ln(s)sds (A.28) 
Jo 

The remaining expressions, Ai,A2,B2 are more complicated and we do not show them 
here. They can be found by expanding cqo in equation flA.lOl) up to fourth order in A, 
as well as coi up to third order in A in eqs.( lA.19l) and ( 1A.18I) . and Cmi up to second order 



for m 7^ and / 7^ 0. Finally, we compute the determinant of the matrix Cmi + ^mi, and 
equate to zero, i.e. 

l + \Bi + \^B2 (A.29) 
+ X\Ao + AAi + X'A^) ln(-exp(7)) = , 

which then directly leads to equation (]A.25p with the identifiable constants Ai, ^2,-82- 
If only Ai is needed, lower orders are required at each step but then ^2,-82 are not 
obtained correctly. 

Appendix B. Matrix Elements 

The different matrix elements with shifted correlated n-dimensional Gaussians can be 
calculated with the help of the following integrals, 

^„^^_,T5,+,T,. ^ expi-v'^B-'v) = M , (B.l) 

V det B 4 

^n^^-x-rsx+v-r^ (fl^x) = {Ju)M , (B.2) 
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{x^Dx) = (^u^Du + ^tT{DB-^)^ M , (B.3) 



^n^^-xTA'x+s'Tx f __f_A_f_ , e 

\ OX ox J 

{2tT{A'AAB-^) + Au^A'AAu - 2u^{A'As + AAs') + s'^As) M , (B.4) 



-x"^ Ax+s^ X 



where 



B = A + A , V = s + s\ u = ]^B . 



(B.5) 



For dipole-dipole interaction equation ([2]) the matrix elements of with the long-range 
Gaussians equation fl2T]) are given as 



e 6^ v/(r, yjjrdrdv? = — — TT 

d 2 



1 



1 rf2 



f/(2,-,-)-2f/(l,-,^,. 



(B.6) 



where 



1 /"°^ 



(B.7) 



is the Tricomi confluent hypergeometric function. Note that only the monopole part, 
Vo(r), gives a non-zero contribution to the matrix element. 
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